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ABSTRACT 

We present a vector formulation of an interferometric observation of a star, including 
the effects of the barycentric motion of the observatory, the proper motions of the 
star, and the reflex motions of the star due to orbiting planets. We use this model 
to empirically determine the magnitude and form of the signal due to a single Earth- 
mass planet orbiting about a sun-mass star. Using bounding values for the known 
components of the model, we perform a series of expansions, comparing the residuals 
to this signal. We demonstrate why commonly used first order linearizations of similar 
measurement models are insufficient for signals of the magnitude of the one due to an 
Earth-mass planet, and present a consistent expansion which is linear in the unknown 
quantities, with residuals multiple orders of magnitude below the Earth-mass planet 
signal. We also discuss numerical issues that can arise when simulating or analyzing 
these measurements. 

Subject headings: astrometry, methods: analytical 



1. Introduction 

Much study has been dedicated in recent years to the possibility of using an ultra-precise, 
space-based interferometer for the purpose of discovering extra-solar planets by their impact on the 
astrometric positions of their parent stars. Th is has become one o f the major science areas of th e 



proposed Space Interferometry Mission (SIM) (jSozzetti et al 



2002 



20031 : ICatanzarite et all 12006) 



and has also been con sidered as an application for the European Space Agency's Gaia mission 
(jCasertano et al.lll996l ). Of particular interest to the exoplanet community is the possibility that 
interferometers capable of sub-//as precision can be used to detect the presence of Earth-sized 
planets in Earth-like orbits — a goal which is many years away from being realized by any of the other 
currently employed or studied planet-finding methods. A number of studies have been completed in 



order to assess the exact planet-fin ding capabilities of astrometric instruments (jTraub et al 



2009 



Casertano et al J 120081 : lBrownll2009i ). One byproduct of these studies has been the r ealizat i on tha t 
the classical description of astrometric observations (as described, for instance, in iGreenl (j 19851 )) 
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makes approximations that are suitable only when dealing with levels of precision of 1 mas or 
higher. Several mor e precise descriptions have been published, including a very thorough one in 



Konacki et al.l (|2002l ) , but most of these take the classical approach of separately treating the effects 



of proper motion, parallax, and the stellar reflex due to companions, with separate expansions of 
each effect. Furthermore, when demonstrating analysis techniques, these studies often still only 
use a first order expansion to simplify the required computations. While it is possible to achieve 
the required numerical precision with these approaches, there is an added burden from having 
to separately consider the expansion of the direction vector and other effects. We believe that a 
simpler approach is to linearize a single measurement equation to produce one unified expression. 

Here, we derive the exactB expression for an astrometric measurement, and then present several 
expansions to multiple levels of precision. This exercise is important for two reasons. First, if one 
wishes to evaluate an algorithm, it is crucial to ensure that any simulated test data does not contain 
biases or components not present in the true data stream. Even if such structures are below 
the level of other simulated noise sources, they may have an effect on any processing algorithm 
which makes the assumption of white, gaussian (or pseudo-gaussian) noise. The added signals will 
not be random, and, as shown below, may closely resemble the signal sought in planet-finding 
applications. For these reasons, we believe that the correct way to simulate astrometric data is to 
use an exact representation of the physical system being modeled. This removes the possibility of 
inadvertently introducing non-random noise sources, or otherwise creating an unfair test for the 
analysis algorithm. Second, when analyzing astrometric data, while linearization of the signal is 
a very useful tool, we must always ensure that such manipulation does not produce a template 
that is measurably different from the data. If data is generated using the same linearization as 
is assumed by the analysis method, and the linearization introduces measurable structure not 
present in the true signal being simulated, then use of the same linearization in both simulation 
and analysis does not constitute a fair test of the algorithm. Therefore, the main focus of this 
paper will not be a specific analysis technique. Rather, we seek to develop an exact formulation of 
the astrometric measurement so that completely unbiased data can be produced, on which various 
analysis techniques can be tested, and to examine the simplifications that can safely be made to 
this exact form. 



*By exact we do not mean that all possible contributions to the measurement are included; for instance, we have 
not yet considered relativistic effects. Rather, we mean that we are formulating the exact form of the nonlinear 
measurement for the given set of effects included: parallax, proper motion, and stellar reflex. 
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2. The astrometric measurement 

Figure Q] and Table Q] define the vectors and reference frames used to describe an astrometric 
observation of an exosystem[j] Reference frame I is an inertially fixed, bary centric/ecliptic frame 
located at the solar system barycenter with the unit vector e3 directed perpendicular to the ecliptic 
plane (the ei and e2 unit vectors are arbitrary). We define a second inertially fixed frame, called 
the tangent frame, B, with b3 axis aligned with f s (io), the unit vector to the star at the initial 
time to- The bi axis is perpendicular to the plane containing r s {to) and 63. The final unit vector 
direction, b2, is mutually perpendicular to bi and r s (to). Given these definitions, we can find the 
ecliptic coordinates in I of the three unit vectors defining B, 



r s (to) = b 3 
bi 
b 2 



cos A cos (3 sin A cos j3 sin (3 

l T 

1 

f s (t ) x bi 



1 T 

1 



x r s (t )/ cos (3 



sin A cos A 



cos A sin (3 — sin A sin (3 cos (3 



(1) 
(2) 
(3) 



where ? s (io) is the unit vector to the star as given by its ecliptic coordinates at epoch to: (A,/3), 
measured in the I frame. This definition of frames assumes constant relative velocities between 
exosystem and local barycenters (i.e., constant proper motions). For the time scales on which 
astrometric observations are taken, this is reasonable, but this derivation may have to be expanded 
if dealing with an accelerating exosystem. 



We define the parallax by the small quantity, 



|r s (*o) 



(4) 



where a is a distance constant such that w is in units of radians. Thus, r s (io) = ( a /w)r s (to) for 
a = 1 AU, with r s (to) expressed in AU. An estimate of w provides the distance to the target star 
at epoch to. For a typical target star at 10 pc (2.0626+ x 10 6 AU)@, ro~5x 10~ 7 . The classical 
astrometric equations retain terms only to first order in w although, as we show in fj3l second order 
terms can be significant when working with //as precision measurements. 

The motion of the barycenter of the target star system, r^, is considered as motion in the 
tangent reference frame, 



= &x(t - to)bi + a y (t - t )b 2 + cr z (t - i )r s (to 



(5) 



^In general, we will use the notation x to denote the unit vector of x (x = x/||x||). All bold symbols refer to vectors 
and terms without explicit time dependence are assumed to refer to time t. We use G to refer to the barycenter of 
the exosolar system (target star and planets). The proper motion is defined as a motion of G. 

'In the simulations in this paper, the conversion between pc and AU is always taken exactly, i.e., 1 pc = 

( tan ( 180x*3600 )) AU - 
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where a% are the components of barycenter velocity at epoch to hi units of distance/time. We 
approximate this velocity to be constant, as is usually done when considering short time spans such 
as a space-based observatory's lifetime. T his expression is generally split into two components: the 



transverse and radial velocities. Following iGreenl (|1985l ). we can write: 

(6) 



gr (t) = S( r G(io) + r,(*)) = 2r M (t) 



V R r s {t ) + V T where V T = r a (t ) X (^r G (t) x f s (i )) 



where the superscript on the derivative denotes differentiation in the inertial frame, and tq is the 
position of the target system barycenter (see Figure [1]). In our notation, Vr is just a z , whereas Vy 
is the vector [a x , a y , 0] T . The transverse velocity (multiplied by time) thus gives the proper motion 
of the barycenter (motion in the plane of the sky) , which has the largest effect on the position of 
the target system at the time of observation. However, the radial velocity also has a measurable 
effect on the astrometric observation. Because motion of the target system in the radial direction 
causes the distance to the system to change, we observe a difference in the direction to the target 
star, known as 'perspective acceleration'. As will be discussed in the third component of r M 
(and thus a z ) interacts in non- negligible ways with other values in the measurement, and is thus 
observable. However, when dealing with a priori estimates of these velocity components, as in £ )3.2l 
it is important to note that the transverse and radial velocities are measured in different ways (i.e., 
astrometry vs. doppler spectroscopy). 

We will find it convenient to have an expression for the normalized barycenter motion, 
r 

f M - Ti — 7TT7T = & x{t ~ *o)bi + a y (t - t )b 2 + a z (t - t )r s (t ) (7) 

where a x , a y , and u z are the nondimensional barycenter velocities in radians/time unit. For typical 
stars considered as candidates for astrometric planet finding, angular velocities from 100 to 1000 
mas are common, making f„, over a 5-10 year period, of order 5 x 10~ 7 to 5 x 10~ 6 . 

Finally, we must consider what exactly is measured during an astrometric observation. An 
interferometer measures the projection of a target direction onto the instrument's baseline vector, 
recorded as the optical path-length delay (OPD) between the two interferometer detectors, 

OPD = B • r s/sc + k + n (8) 

where B is the orientation vector of an interferometer of baseline length B = ||B||, k is a constant 
term representing the offset of the optical path differences, and n is the measurement noise. Planet- 
finding is generally proposed in a narrow angle mode, where the measurement is the relative OPD 
between two sources in a field of view, performed in quick succession such that k remains nearly 
constant, and is assumed to cancel in subtraction. Typically each target star has multiple reference 
sources, so the combined differential OPD (AOPD) is with respect to a centroid position. The 
measurement is also repeated for two (preferably orthogonal) interferometer baseline orientations 
to track the 2D position of the target in the plane of the sky. Assuming that the interferometer 
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baselines are taken as directions t>i and i>2 of our body frame, then the astrometric measurement 
becomes: 

b 2 • fts/sc - r c / sc ) 

where B is the size of the interferometer baseline and n is an additive noise vector due to measure- 
ment error. 



B 



+ n 



0) 



In the literature, the measurement in equation ([9]) is often described as the angular separation 



betw een the two sources, projected onto the baseline (jKonacki et al 



2002 



Sozzetti et al. 



2002, 



2003); it is assumed that d scaled by B is a radian measure that can be converted to other 
angular units such as arcseconds. It is this assumption that leads to the terminology of /uas-precise 
astrometry, as the required sensitivity of the instrument to changes in the OPDs, when treated as 
an angle, evaluates to under 1 /xas. 



In fact, IColavital (|1994l ) points out that, using the definition for r s / sc in equation for a 
centroid separation of (AA, A/3), the difference in unit vectors, to first order, can be written as, 



^s/sc '-c/sc 



sin j3 cos AA/3 + cos /3 sin AAA 
sin f3 sin AA/3 — cos (3 cos AAA 
-cos/3A/3 



(10) 



In this way, knowledge of the baseline vector and the differential OPD allows you to calculate a 
vector which maps in a relatively simple fashion (to first order) to two spherical angles representing 
the separation between the target and centroid. 

Unfortunately, such a direct first order mapping between equation ([9]) and the angle between 
the target and centroid is not very accurate. An OPD has units of distance or time (with the 
converting factor equal to the speed of light), and so a differential OPD will also have values equal 
to fractions of the interferometer baseline. Let us assume that the (flat) wavefront from the target 
star is incident on the interferometer baseline with angle 6i, and that the centroid is located at an 
angle A9i from the target star, in the projection of the sky due to interferometer orientation 
(Figure[2] — A6i is analogous to AA and A/3 in equation (|10p ). The components of the measurement 
(without the noise term) are then, 

di = B (cos 6i- cos(6i- A6i)) = B (cos 6i(l- cos A6i) -sin 8i sin A8i) . (11) 



We are primarily interested in changes in the angle A6i, since this determines the movement 
of the target star with respect to the centroid in time (of course, this is further complicated by 
possible motion of the centroid itself). Assuming this to be a small angle, we can substitute the 
Taylor series expansions of the sine and cosine terms in A6i to first order to find, 

di » -BsinOiAOi. (12) 

Thus, by scaling the differential OPD by the interferometer baseline length and the direction 
on the sky of the target (-Bsin^j), we do get a first order approximation of the angular difference 
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between the target and centroid (assuming perfect a priori knowledge of the target's (or equivalently 
centroid's) location with respect to the interferometer ori entation). Unfortunately, if we assume 
the t arget star to be on the order of 1° from the centroid (jSozzetti et al.l 12002 ; IShao and Colavita 
19921 ). the first of the dropped terms (AO?/ 2) has a magnitude of 1.5xl0 -4 rad, or 31.4 arcseconds. 
When discussing ultra-precise applications, it is therefore inaccurate to treat the values produced by 
equation ([9]) as directly mapping to angular measures. For these reasons, we treat all derived values 
as dimensionless for the remainder of this discussion, normalizing all distances and converting all 
angles to radians to remain consistent. It is also important to point out that in this discussion, we 
consistently assume zero pointing error. For a real instrument, accurate knowledge of the baseline 
vectors will be built up over many individual observations, each with an associated measurement 
error, but neither the interferometer orientation, nor the exact length of the baseline, can ever be 
known with perfect precision. In order to include a pointing and baseline length error, we would 
need to extend the measurement equation, and update all subsequent calculations. 

With these definitions and considerations, we can write the exact astrometric measurement 
in terms of known quantities and the quantities whose values we wish to determine. This is done 
by writing the vector from the spacecraft to the target star, r s / sc , in terms of the initial reference 
vector, r s (t ), 



?s/sc = r s (t ) - r s/G (t ) + +r s/G - r s 
= r s (t ) + r^+Ar s/G -r sc . 



(13) 



where r sc is the spacecraft position vector relative to the solar system barycenter and Ar s / G is the 
difference in the star's position relative to G between to and epoch. We can then find the unit 
vector in the direction of r s / sc , 



s/sc 



s/sc 



W^s/scW 

(r s (t ) + r /i + Ar s/G - r sc ) x 

r a (i ) • r a (t ) + !> • r M + Ar s/G 



+ 2r s (t ) t„ 



+2r s (t ) ■ Ar s/G - 2r s (t ) ■ r sc + 2r^ • Ar s/G - 2v fl • r sc - 2Ar s/G ■ r sc 



.(14) 



The spacecraft to centroid pointing unit vector, r c / gc , can be similarly expressed by using 
equation f|X4j) for each reference source with the appropriate values for r s (to) and r„. When using 
this model for the reference sources, it is important to note that for a source n, f n (to) is not 
aligned with r s (to) and is thus not orthogonal to bi and b2. When evaluating equation (jUJ), this 
simply means that all vectors have to be expressed as components in the same reference frame, 
which requires a change of coordinates for the reference sources. This does, however, lead to some 
complications when dealing with linearizations of the model, which will be discussed further in 
Jj3l If we assume extragalactic references with no companions or planets (which was done for the 
simulations described in this paper) the variation in the centroid position will be many orders 



-7- 



of magnitudes below our desired level of accuracy. Unfortunately, such references are not always 
available, and so these effects may become significant if one or more of the references is within our 
galaxy and has companions of its own. 

When simulating an astrometric data set, it seems natural to use the formulation in equation 
(I14p . The effects of parallax and proper motion are represented exactly, with no ambiguities in 
how perspective acceleration should be introduced. Furthermore, the coupling of these effects with 
the astrometric wobble due to planetary systems also appears and does not need to be separately 
considered. The computational effort required by this equation is not significantly greater than 
that of the various measurement expansions commonly in use, and is actually less than that of the 
second order expansions presented in the next section. Other than a few numerical considerations 
discussed in §41 the application of this formulation is trivial, and ensures high confidence in the 
simulated data. 



3. Measurement Expansion 

While equation (j!4h provides a method for exactly simulating precise astrometric measure- 
ments, the analysis of such signals can be made difficult by the nonlinear interactions of the various 
terms in the normalization factor. The classical approach has been to expand this equation and 
retain only relevant terms. As we will show below, for the purposes of planet-finding, these must 
necessarily include nonlinear terms, but an expansion is still useful as a simplifying step in the data 
analysis. In order to decide which terms should be retained, we must first quantify the desired 
precision of the measurement. There are several ways to do this, but a relatively straightforward 
one is to calculate the magnitude of the smallest signal we wish to be able to measure. As our 
focus here is on the detection of exoplanets, we will take as our smallest signal the effect of a single 
Earth-mass planet on an Earth-like orbit on the astrometric signature of a sun- twin star. Figure [3] 
shows the differences between the results of two simulations, each using equation (I14p and a fixed 
centroid position (with respect to the solar system barycenter). The simulations include identical 
values for parallax, barycenter motion, and assume a sun-twin target star lying exactly 10 parsecs 
from the solar system barycenter (at epoch to). One simulation, however, includes the stellar jitter 
due to an Earth-mass planet, while the other does not (the line of nodes of the planet's orbit is 
rotated by 45° with respect to the line of sight). The difference between the two measurements is 
the most basic representation of the magnitude of the signal in which we are interested, and, in 
this case, shows a signal on the order of lxlO 12 . This value is in the units of equation ([9]): it is the 
difference between fractions of the baseline distance. 

This tells us that to be sensitive to the jitter due to an Earth- twin with any measure of 
confidence, we must retain all terms in the measurement of order 1 x 10 13 or greater (assuming that 
our analysis technique will be sensitive to structured signals up to one order of magnitude below 
the target signal). To aid in this calculation, Table [2] enumerates the ranges of typical values for 
the various terms in equation (114p . 
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3.1. Expansion Assuming No Prior Knowledge 

We begin the expansion by dividing the numerator and denominator of equation (|14p by 
|r s (to)|| to find the normalized version, 

h/sc = (r s (*o) + f M + roAf s/G - wr sc ) x (15) 

1 + • f M + tu 2 Af s/G • Af s/G + tu 2 f sc • r sc + 2? a (t Q ) • 
+2wr s (t ) ■ Ar s/G - 2wr s (t ) ■ r sc + 2ror Ai • Af s/G - 2wf^ ■ r sc - 2w 2 Ar s/G ■ f sc 



where w = a/||r s (io)|| and r^ = r M /||r s (io) ||, ? sc is the spacecraft position normalized by a and 
Af s / G is Ar s / G , also normalized by a. 

We next apply a binomial expansion to the denominator of equation (]15p and retain terms 
explicitly in w to first order. Since the signal we are looking for appears explicitly only as GjAf s / G , 
any terms of proportional order must be retained. Thus, terms proportional to ||f M || 2 must be 
left in, but terms proportional to Ijf^ll" for n > 2 and terms proportional to Hf^l^tu for n > 1 
are dropped. We also drop all terms proportional to r s (to), as dotting r s / sc with bj as in the 
measurement equation causes these terms to equal zero. The resulting approximation is: 

f • b- w [ ifl + tuAfs / G + ro (^(*o) • f sc )f M - ro(Af s/G • r s (t ))r M - wi sc \ fe 
S/SC 1 \ - (r s (t ) -f M )f M + ti7(r s (to) •f M )f sc -ro(r s (t ) •f A1 )Af s/G J 

This is a linear expression in parallax and can be used to extract the stellar reflex signal, though 
it contains complicated couplings among the spacecraft position, parallax, barycenter motion, and 
stellar reflex. Note that it is not assumed that these other terms are known; rather, the barycenter 
motion and the parallax factor, w, must also be estimated in the data analysis. This leads to the 
problem of attempting to fit the radial motion of the target from purely astrometric measurements. 
While the radial and proper motions of targets are approximated as completely separate in the 
classical treatment, the level of precision required for this application means that we must consider 
the effects of radial motion. The first term of equation [TB] is f„ • bj. If bj is, as we have assumed, 
equivalent to bi or b2, then the a z dependence of this term is exactly zero. Similarly, the next two 
terms in the direction of f« will have no a z components when dotted with bj. However, the two final 
terms in this expansion have magnitudes given by va(r s (to) which makes them proportional to 
a z . As these two terms are in the directions of f sc and Af s / G , which are arbitrarily oriented in the 
tangent frame, this a z dependence is not zeroed even when the baseline is perfectly aligned with 
bi or b2. Thus, the radial star motion explicitly enters our expression, and makes a significant 
contribution at the desired level of precision. 

Figure H] repeats the simulation from Figure [31 including the planet in both cases, but now 
comparing the exact expression for r s / sc and the first order expansion in equation [TBI The resulting 
error is less than one order of magnitude below the signal due to the planet, which means that it is 
of the same order as the desired sensitivity for all targets closer than 10 pc or with any combination 
of the other parameters which would cause the planet signal to decrease in magnitude. Additionally, 
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the structure of the error is dominated by parallax terms, which produce periodic structure that 
could be mistaken for astrometric wobble due to another planet. We can compare equa tion [TBI 



with previously published first-order forms such as Equations 31 and 32 in lKonacki et al.l (|2002l ). 
These have a very similar form to equation [TB| save that it is assumed that all radial components 
of motion are unobservable, and thus the equivalent vectors to Ar s / G and are formed without 
the r s (to) components. The expression there also does not subtract the initial displacement of the 
target star from its system's barycenter, thereby applying the proper motion to the star rather 
than the system barycenter. If we redo our simulation using the first order expansion with the 
radial terms dropped, we end up with residuals one order of magnitude higher than those in figure 
[H Thus, the exclusion of radial terms produces residuals that are one order of magnitude above 
the target signal, rather than one order of magnitude below, as seen in Figure [5j 

Furthermore, a number of the terms dropped in this expansion were proportional to zu 2 multi- 
plied by quantities of order 1. The astrome tric literatur e is quite clear that a first order expansion 



in zu is only good to milliarcsec accuracy. (|Greenlll985l ) If our instrument has a final precision of 



sub-/iarcsec, it is important to ask if measurable terms were dropped in the expansion. Such terms 
can get as large as 6 times zu 2 . Again, for a star at 10 pc, ro~5x 10~ 7 , making the second order 
terms of order 2.5 x 10 -13 . Thus, for most stars the approximation is a good one but it does raise 
concerns for the closest targets, or for analysis methods that are sensitive to non-random signals 
just below the level of the target signal. 



To address this, we can repeat the expansion to second order in zu. However, for consistency, 
we must now retain some additional terms. The zeroeth and first order terms (in zu) remain the 
same (since 1 1 r ^ 1 1 3 is of order lxl0~ 16 and ||f^|| 2 ro is of order lxlO -17 ). Any terms in zu 2 with 
factors of can safely be dropped. Terms in zu 2 with factors of Af s / G should also be quite small, 
but we will leave these to first order since it is theoretically possible that there exist systems with 
both Earth-like planets and very large, widely separated super-Jupiters which have not yet been 
detected by RV (due to their very long periods), and which would cause a large stellar reflex to 
make these factors significant. These considerations leave only three new terms in the expansion: 



/ r M -(r s (i ) 



■b, 



zu 



\ 



+ ((? s (i ) • f s 
+ (r sc (r s (t ) 



+ (Ar s/G - r sc - (r s (t 
~ ( Af s/G • f s (t ))f^ + (r s (t ) ' f M )f sc ) zu 



r M )Af s/G ) 



•b; 



(17) 



Ar s / G ) + Af s / G (? s (i ) • f sc ) - r sc (r s (i ) ■ ? sc )) ^ 2 



Figure E] shows the difference between the exact expression for r s / sc and the second order expansion 
in equation fllTf) . The resulting error is now two to three orders of magnitude below the target signal. 
The error structure is now dominated by proper motion terms, which are linear, not periodic. These 
could still be confused with the effects of long period planets, but the signals are so low compared to 
the target signal that this possibility is unlikely. It is important to note, however, the fitting problem 
has become significantly more challenging as the parameters enter nonlinearly in the measurement. 



As has already been mentioned, our practice of dropping terms that are in the direction of 
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?s(io) breaks down in two cases: when there is pointing error such that the interferometer baseline 
is not aligned with bi or l>2, or when expressing the direction to a different star (i.e., a non- ideal 
reference source) in the tangent frame defined by the target star. In these number of 

additional terms must be included in our expansion. Following the same procedure as above, we 
can rewrite equation (fT7|) as: 



' r a (t ) + r fl - t s (t )(r a (t ) ■ 17, 



~s/sc 



-r M (r s (i ) t^) 



\ 



+ 



V 



&r s/G - r sc - Af s / G (r s (t ) ■ r M ) + ? sc (r s (i ) ■ r M ) 
- f s (i )(f s (i ) • Ar s/G ) - f M (r a (t ) ■ Af s/G ) 
+ 3f s (t )(r a (t ) • f M )(r s (i ) • Ar s/G ) + f s (t )(r s (t ) ■ ? sc ) 
(r s (i ) • f sc ) - 3f s (io)(f s (io) ■ f^)(r s (t ) ■ f sc ) 



+ r 



- r s (t )(f M • Ar s/G ) + r s (t o )0 
r sc (r s (t ) ■ Ar s/G ) + Af s/G (f s (t ) • f sc ) - r sc (r s (t ) ■ f sc ) 

- 3f s (t )(f s (to) • Af s / G )(f s (t ) • f sc ) + f s (t )(Af s/G • f sc 



ZD 



ZD' 



(18) 



where the subscript s may now refer to either the target star, or any of its reference sources. 

We now have a completely general second order expansion of the spacecraft to source unit 
vector, which applies when the source is different from the target star used to define the tangent 
frame, or when the spacecraft baseline pointing does not coincide with the vectors of the tangent 
frame. In the first case, we retain terms proportional to f s (to) ' hi because r s (to) 7^ b3 and so 
these terms do not automatically go to zero. In the second case, these terms do not go to zero 
because bj 7^ bi or b2. This expression is even more complex than equation (|17p . and even less 
suitable for use in analysis. However, this added complexity is actually mitigated by the same 
reasons that caused us to add in the extra terms. First, if the pointing error is small (i.e., the 
baseline is reasonably close to one of the two tangent frame unit vectors) , then only the new terms 
in r s (to) whose coefficients are of order 1 should be kept, as f s (to) ' b« will still represent a small 
value. Second, even non-ideal reference sources will still usually be stars which are significantly 
further from the solar system than the target star, with much smaller parallax and barycenter 
motion values. As such, the first order form will typically be sufficiently precise for these sources, 
allowing most or all of the terms in zu 2 to be dropped. 



3.2. Expansion With Prior Knowledge 

While the expansions in £ j3.1l are a useful representation of an astrometric measurement, they do 
not quite correspond to the real-life analysis problem of astrometric planet-finding. In fact, before 
any high precision astrometric measurements are made, we will already have some knowledge of 



-li- 



the barycenter motion and parallax of our target stars, which will only need to be updated by some 
small amount. To represent this, we can rewrite equation (|15p with f « replaced by r^o + 5fn and m 
replaced by vjq + 5m, where r^o and vjq are the currently known components of barycenter motion 
and parallax and 5r„ and 5m are the (relatively small) correction terms. The result is: 

r s /sc = (f s (io) + *>o + + ( w o + <^)Ar s/G - (w + 5m)r sc ) x (19) 

— - 

1 + f/io • f M o + ^ " + 2f M) • Sf,j, + (m + <5to) 2 Ar s/G • Af s/G + (ro + ^ro) 2 f sc • f sc 
+ 2r s (i ) ■ f M o + 2f s (i ) • Sf^ + 2(tu + 5m)i s (t ) ■ Af s / G - 2(ro + <ku)r s (t ) ■ f sc 
+ 2(zu + <5ro)(r M o ■ Ar s/G + ^ • Af s/G ) - 2(w + <5ro)(f At0 • r sc + Sf^ ■ r sc ) 
- 2(zu + 5m) 2 Af s/G • f sc 



In order to carry out expansions as before, it is necessary to quantify the expected values of 
these new variables so that we can decide which terms are negligible and which are important. 
Since we are assuming that our a priori estimates will be reasonably good, r^o and mo should be 
of the same order as and m above. Current ly, parallax and pr oper motion are known to at least 



1 mas and 1 mas/year accuracy, respectively (jTuron et al.lll995l ). We will take 5m to be of order 



5x10 9 ; ||#r„|| will vary in magnitude in time, but could be of order 1x10 8 within five years of 



observations. 

We again perform a binomial expansion on the denominator of equation (|19p and retain terms 
explicitly in wq to second order and explicitly in 5m to first order. As before, terms proportional 
to Hf^oH" for n > 2 and terms proportional to ||r„o|| n ^o f° r n > 1 are dropped, as are terms 
proportional to ||<5f M || n for n > 1, terms proportional to ||<5r M ||<5ri7, and terms including m 2 that are 
proportional to ||A? S / G || n for n > 1. 

Terms proportional to ||f M o||^, ||f || II^uIIj an d ||<5f M j|tz7o are borderline — for the values as- 
sumed here, they have magnitudes of up to lxlO -14 , which means they should be neglected, 
but if the proper motions are larger than expected (or there is higher uncertainty in the a priori 
measurements), then these terms can become significant, so we will leave them. Similarly, terms 
proportional to m<^5m have magnitudes of lxl0~ 15 for the values assumed here, and are fairly 
unlikely to become significant unless uncertainties in parallax are much greater than expected for 
the nearest stars. We include only one of these terms because leaving it out leaves a clear sinusoidal 
pattern in the signal of order 5x 10~ 15 . If the chosen analysis method is not expected to be sensitive 
to such a signal, this term can safely be omitted. The resulting approximation is: 
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'■"s/sc 



( r M o + 5r^ - r M o(r s (to) • r M o) - <h>(r a (i ) • r M0 ) - r M o(r s (t ) ■ <5r M ) ^ 
Ar s/G -f sc - A? s/G? (f s (£o) -f M o) 1 , 
+ r sc (r s (t ) • i>o) + r> (r s (io) ■ r sc ) _ 

- r^oCfs^o) • Af s/G )(5ro + [f sc (r s (i ) • 5f^) + <5r M (r s (t ) ■ f sc )] ro 

- 2f sc (f s (t ) • r sc )ro 5ro 
\ + [r sc (r s (t ) • Ar s/G ) + Ar s/G (i s (t ) ■ r sc ) - r sc (r s (t ) ■ r sc )] w\ J 



The difference between this equation and the OPD using the exact expression for f s / sc is 
almost exactly the same as for the second order expansion from £13.1i In fact, the difference between 
differential OPDs calculated using equations (fTT)) and (f20|) is actually at the level of the numerical 
precision of the data type used in these simulations (Figure [7]). However, unlike equation (|17p . this 
equation is linear in the quantities to be estimated, which can make analysis much simpler. 

One assumption made in the above analysis is that the error is of the same magnitude for all 
three components of r„. However, this ignores the fact that currently existing estimates for the 
first two components (the transverse motion) are derived from different sources than estimates for 
the third component (the radial motion). While proper motions have been measured by space- 
based astrometric instruments, radial velocities are generally provided by ground-based doppler 
spectroscopy, and are inherently less accurate. To address this, we can modify equation ([20|) by 
separating 5r^ into two separate errors so that: 

*n = ffjo + 5r R + 5r T (21) 

where 5r R and 5r? are due to errors in our prior knowledge of Vr and Vr, as defined in equation 

We can now repeat the expansion performed to generate equation (|2"U|) . substituting the right- 
hand side of equation f)21 1) for r„. The rules for keeping terms containing 5tt remain the same as 
those for 5f„. If we assume that <5r# is less than three orders of magnitude higher than Sr? (i.e., m/s 
radial velocity precision), then the expansion in equation (|20p must also retain terms proportional 
to ||<5rR|| n for n <= 2. The result is that only one term, equal to — 5rji(r s (to) • Str) = —Srji\\Srji\\, 
must be added to equation (|20p to account for the lower precision of the prior radial velocity 
measurement. 



4. Numerical Considerations 



The final simulation results presented in Figure [7] highlight a very important point. Due 
to the small magnitudes of many of the quantities utilized in this analysis, when it comes time to 
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numerically manipulate them, the limitations of the numerical representation become a real concern. 
Any computer generated value will have an associated finite precision, but specific representation 
schemas may actually introduce precision differences as functions of value magnitude. For example, 
in MATLAB (which was used for the simulations in this paper), the default numerical data type is 
a double precision, floating point number, encoded as, 

2 e x s 

where e is the exponent and s is the significand. These are stored as 11 bits and 53 bits, respectively, 
as per the IEEE 754 standard. Because the exponent and significand are encoded with a finite 
number of bits, the minimum step size to the next available value in this schema is a function of 
the curren t value . Thu s, differences between pairs of values of different orders will have different 



precisions (|Molerll2004l ). For example, when encoding values of order 1 as double precision floats, 
the minimum step size to the next available value is on the order of lxlO -16 . However, when 
encoding values of order le6 (like, for example, the value of 10 pc in AU), the step size to the next 



available value is of order 1x10 10 . 

This becomes important when one wants to generate 'ground truth' data, using, for example, 
the exact representation of the astrometric observation in equation (|14p . In this case, certain values, 
such as v s (to) are multiple orders of magnitude greater than 1, whereas the final result is many 
orders of magnitude less. This means that numerical noise at the level of the desired signal could 
inadvertently be introduced into the simulated data. In the simulations presented here, this only 
became significant when converting r s j sc to a unit vector. One solution (the one employed here) is 
to generate the large values with an arbitrary precision data type using the appropriate number of 
digits. Once the spacecraft to star unit vector is generated, all values are of sufficiently low order 
that the default fixed precision data type is sufficient for the remaining calculations. To do this, it is 
necessary to use an arbitrary-precision data type, which simply encodes values using variable length 
significands, rather than the fixed bit scheme described above. For the simulations here, this was 
done using the GNU multiple precision library (http://gmplib.org/), with the significand stored 
as 256 bits. The numerical effect is quite limited in the simulations above (with the target star at 
10 pc), but becomes much more pronounced for targets at 100 pc. Figure [8] shows the difference 
between differential OPDs using equation (|14p . with one simulation done using the default double 
precision data type for all values, and the other using the multiple precision data type. The error 
between the two is close to 1 x 10~ 8 in magnitude, making it a significant noise contribution. 

An alternative to this approach would be to use equation (|15|) . All terms in this formulation 
(which is equivalent to equation (fT4"l) ) are of order 1 or less, so that the default double precision data 
type should be sufficient. However, use of this equation can also introduce numerical noise, if one 
is not careful in how the various terms are defined. Terms such as are essentially small values 
divided by very large values. If they are not exactly defined, but rather calculated in software, 
then they will have the same problems as shown above. Since this discussion applies only to the 
simulation or generated data, it is relatively simple to just assume values for barycenter motions 
that will give exact fractional values when scaled by the parallax, eliminating such concerns. 
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Finally, if we wished to update our expansions from $3] such that the residual was below the 
precision of the default double precision data type, we would need to retain terms proportional to 
tu 3 . However, in the case of the expansions assuming prior knowledge of parallax and barycenter 
motion, no higher order terms in 5m would be needed. 

5. Conclusions 

In this treatment, we have presented an exact vector formulation of a narrow-angle astrometric 
observation incorporating the effects of parallax, stellar reflex, and barycenter motion. Because the 
derivation leading to equations (fT4"j) and (fT5|) does not represent a significant computational burden 
for modern computers, we believe that they (or something fairly similar) should be used for the 
simulation of data for the purpose of testing analysis techniques, rather than any linearization of 
any order of precision. While it may be argued that using a simplified data set (i.e., one derived 
from approximate or linearized descriptions of the true observation) allows for initial testing of an 
analysis method, which can then be followed up with more refined tests, we do not believe that such 
experiments are scaleable. The most important point to keep in mind when using astrometry for 
planet finding is that the expected signals will be of very low order, and will interact non-linearly 
with several unknown (or partially known) parameters. Therefore, testing data produced by a 
linearization which may itself introduce residual signals close to the order of the planet signal, and 
which adds the planet signal without modeling its interaction with other terms, does not actually 
tell us anything about our analysis technique's ability to isolate the planet signal in a real data set. 

Another argument which is often put forward is that such considerations are unimportant 
because the measurement noise will generally be higher than the errors discussed here. However, 
most analysis techniques (be they Bayesian inference or nonlinear programming minimization), 
make some assumptions as to the structure of the noise. Even if the noise is not described as additive 
white Gaussian, any autocorellation and non-zero mean components are carefully modeled based on 
our knowledge of the physics of the observed systems and instrument. Linearization residuals are 
not random and often introduce patterns that can be quite similar to planet signatures, especially 
due to the parallax effect of any sun-orbiting observatory. Use of the exact observation description 
obviates the need for all such considerations. 

Furthermore, the linearization presented in §3.21 allows us to use a sufficiently precise linear 
expression for analysis, with only the assumption that we have measurements of parallax and 
barycenter motion of certain fidelity. This should be very helpful for any methods reliant on 
fitting techniques, as the nonlinearities in the second order expansion from i j3.ll are quite difficult 
for most least-squares algorithms. Again, both the functional representation and computational 
requirement of this linearization are only slightly more demanding than those of the classical first 
order astrometric equations, and so there appears to be no reason not to use representations of the 
astrometric observation derived in a way similar to the one shown here. 
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Table 1: The vectors used to determine the astrometric measurement. 
Symbol Definition 

r G(to) Vector from solar system barycenter to G at epoch to 

r sc Vector from solar system barycenter to spacecraft at time t 

r pm Vector from solar system barycenter to G at time t (movement due to barycenter motion) 
ippm Vector from spacecraft to G at time t (movement due to barycenter motion and parallax) 

Vector from G(to) to G at time t (barycenter motion) 
y s /q Vector from G to the star at time t (equal to zero if star has no planets) 
r s Vector from solar system barycenter to the star at time t 

r C/ / sc Vector from spacecraft to centroid of reference stars 
r s / sc Vector from spacecraft to star at time t (the measurement vector) 



Table 2: Typical values for terms in equation (|14p . 



Term 


Typical Order 


Notes 


a 


1 


All distances are assumed to be in AU 


r a (*o) 


1 




— * 


5xl0 _6 /year 


Barnard's star is 5xl0~ 5 /year 




lxKT 6 to lxl(T 2 


For solar system planets 


'"sc 


1 


Assuming a ground-based, Earth orbiting or Earth-trailing spacecraft 


W 


5xl0~ 6 to 5xl(T 8 


For target distance of 1-100 parsecs. 



* While r M is a time dependent value, for measurements made over a span of 5 to 10 years it will change by less than 
one order, so we will assume a conservative, constant value when considering which terms to keep in our expansions. 



This preprint was prepared with the AAS L^TgX macros v5.2. 
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Fig. 1. — A schematic of an astrometric observation. The solar system barycenter is placed at 
the origin of frame X, while the observed star system's barycenter at epoch to is at G(to) and at 
G(t) at observation time t. Point s represents the star position at time t and point c represents the 
(possibly time- varying) position of the centroid of a group of reference stars. 




Fig. 2. — A schematic of a narrow- angle astrometric measurement. B is the size of the interfer- 
ometer baseline. 
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Fig. 3. — Difference between differential OPDs (between a sun-twin star and fixed centroid) using 
exact expressions for r s / sc with and without the effects of an Earth-mass planet on an Earth-like 
orbit. The two curves represent the measurements along two orthogonal interferometer baseline 
orientations. 
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Fig. 4. — Difference between differential OPDs (between a sun-twin star and fixed centroid) using 
the exact expression for r s / sc and the first order expansion from equation (|16f) . The two curves 
represent the measurements along two orthogonal interferometer baseline orientations. 
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Fig. 5. — Difference between differential OPDs (between a sun-twin star and fixed centroid) using 
the exact expression for r s / sc and the first order expansion from equation (|16p with the radial 
terms neglected. The two curves represent the measurements along two orthogonal interferometer 
baseline orientations. 
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Fig. 6. — Difference between differential OPDs (between a sun-twin star and fixed centroid) using 
the exact expression for r s / sc and the second order expansion from equation (|17p . The two curves 
represent the measurements along two orthogonal interferometer baseline orientations. 
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Fig. 7. — Difference between differential OPDs (between a sun-twin star and fixed centroid) using 
the second order expansion of i s / sc from equation (|17p and the expansion from equation (|20p with 
assumed errors in parallax and barycenter motion of 1 mas and 1 mas/year, respectively. The 
pattern in the residuals is due to the precision of the numerical data type used. The two curves 
represent the measurements along two orthogonal interferometer baseline orientations. 
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Fig. 8. — Difference between differential OPDs (between a sun-twin star and fixed centroid) using 
the exact formulation of r s / sc from equation (fl4j) with the default double-precision data type and a 
multiple precision data type at 256 bits. All values are the same as in previous simulations, except 
that the target star is placed at 100 pc. The two curves represent the measurements along two 
orthogonal interferometer baseline orientations. 



